import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
import sys
from sklearn.metrics import roc_auc_score

sys.stdout = open("/REDACTED/NC_BIFSG_metrics.txt", "w")

## Linear Estimator
def regEstimate(dataset, pbvarb = 'predicted_prob_black', outcome = 'audited', wvarb = None):
	if wvarb is not None: 
		model =  smf.wls(outcome + ' ~ ' + pbvarb, dataset, weights=dataset[wvarb]).fit(cov_type='HC1')
		coef= model.params[pbvarb]
		se = model.bse[pbvarb]
		return coef, se, model
	else:
		model =  smf.ols(outcome + ' ~ ' + pbvarb, dataset).fit(cov_type='HC1')
		coef= model.params[pbvarb]
		se = model.bse[pbvarb]
		return coef, se, model

### Operational audits and predicted race probabilitiess
dataBISG = pd.read_csv("/REDACTED/data/final/individualBISG2014_full_final.csv")
dataBISG= dataBISG[~dataBISG.predicted_prob_black.isna()]

nc_BISG = dataBISG[dataBISG['state']=="NC"]
nc_BISG['pprob_black'] = nc_BISG['predicted_prob_black']
nc_BISG = nc_BISG[['taxpayer_id', 'pprob_black']]

### NC Data
##NC data
ncdata = pd.read_csv("/REDACTED/NC_Analysis_Dataset_2014_tplevel.csv")

nc_weight = pd.read_stata('/REDACTED/BIFSG/nc_rewt_2014_coded_output_v4_December_ncsamp.dta')


##########################################################
#### Table A8
##########################################################

ncdata = ncdata[ncdata.black_ind.notna()]
nc_weight = nc_weight[['taxpayer_id', 'black_prob', 'uswgt']]

ncdata = pd.merge(ncdata,nc_weight ,on='taxpayer_id')

ncdata = pd.merge(ncdata,nc_BISG,on='taxpayer_id')
ncdata.drop(['predicted_prob_black', 'black_prob'], axis = 1, inplace = True)
ncdata['predicted_prob_black'] = ncdata['pprob_black'] 
ncdata.drop(['pprob_black'], axis = 1, inplace = True)
ncdata['p_black_rd'] = ncdata['predicted_prob_black'].round(2)

ncdata['unitwt']=1
ncdata['aud_no_research_audits'] = [1 if (x.find('[80]') == -1)
                         and (x.find(' 80]') == -1)
                         and (x.find('[80 ') == -1)
                         and (x.find('[91]') == -1)
                         and (x.find(' 91]') == -1)
                         and (x.find('[91 ') == -1)
                         and y == 1
                         else 0
                         for x, y in zip(ncdata.audit_source_code.astype(str), ncdata.audited)]

ncdata['aud_no_research_audits']=pd.to_numeric(ncdata['aud_no_research_audits'])
ncdata['predicted_prob_black']=pd.to_numeric(ncdata['predicted_prob_black'])
ncdata['p_black_rd']=pd.to_numeric(ncdata['p_black_rd'])
ncdata['uswgt']=pd.to_numeric(ncdata['uswgt'])
ncdata['unitwt']=pd.to_numeric(ncdata['unitwt'])
ncdata['black_ind']=pd.to_numeric(ncdata['black_ind'])


ncdata['pb_50'] = np.where(ncdata['predicted_prob_black'] >= 0.50, 1, 0)
ncdata['pb_75'] = np.where(ncdata['predicted_prob_black'] >= 0.75, 1, 0)
ncdata['pb_90'] = np.where(ncdata['predicted_prob_black'] >= 0.90, 1, 0)

ncdata['pnb_50'] = 1 - ncdata['pb_50']
ncdata['pnb_75'] = 1 - ncdata['pb_75']
ncdata['pnb_90'] = 1 - ncdata['pb_90']


### EITC pop
eitc_ncdata = ncdata[ncdata["eic_ind"]==1]


## recalibrate prob_black
B_b_model = sm.OLS.from_formula("black_ind ~ predicted_prob_black", data = ncdata).fit(cov_type = "HC1")
ncdata['bhat'] = B_b_model.predict(ncdata['predicted_prob_black'])

ncdata['pbhat_50'] = np.where(ncdata['bhat'] >= 0.50, 1, 0)
ncdata['pbhat_75'] = np.where(ncdata['bhat'] >= 0.75, 1, 0)
ncdata['pbhat_90'] = np.where(ncdata['bhat'] >= 0.90, 1, 0)

ncdata['pnbhat_50'] = 1 - ncdata['pbhat_50']
ncdata['pnbhat_75'] = 1 - ncdata['pbhat_75']
ncdata['pnbhat_90'] = 1 - ncdata['pbhat_90']

## recalibrate for eitc
B_b_model_eitc = sm.OLS.from_formula("black_ind ~ predicted_prob_black", data = eitc_ncdata).fit(cov_type = "HC1")
eitc_ncdata['bhat'] = B_b_model_eitc.predict(eitc_ncdata['predicted_prob_black'])

eitc_ncdata['pbhat_50'] = np.where(eitc_ncdata['bhat'] >= 0.50, 1, 0)
eitc_ncdata['pbhat_75'] = np.where(eitc_ncdata['bhat'] >= 0.75, 1, 0)
eitc_ncdata['pbhat_90'] = np.where(eitc_ncdata['bhat'] >= 0.90, 1, 0)

eitc_ncdata['pnbhat_50'] = 1 - eitc_ncdata['pbhat_50']
eitc_ncdata['pnbhat_75'] = 1 - eitc_ncdata['pbhat_75']
eitc_ncdata['pnbhat_90'] = 1 - eitc_ncdata['pbhat_90']

print("Full Population\n")

print("Area under ROC curve")
roc_auc_score(ncdata['black_ind'], ncdata['predicted_prob_black'])

print("\n50% Threshold")
print("False Positive")
ncdata[ncdata.black_ind == 0].pb_50.mean()
print("True Positive")
ncdata[ncdata.black_ind == 1].pb_50.mean()
print("False Negative")
ncdata[ncdata.black_ind == 1].pnb_50.mean()
print("True Negative")
ncdata[ncdata.black_ind == 0].pnb_50.mean()
print("Precision")
ncdata[ncdata.pb_50 == 1].black_ind.mean()
print("Recall")
ncdata[ncdata.black_ind == 1].pb_50.mean()
print("Accuracy")
len(ncdata[ncdata.black_ind == ncdata.pb_50])/len(ncdata)

print("\n75% Threshold")
print("False Positive")
ncdata[ncdata.black_ind == 0].pb_75.mean()
print("True Positive")
ncdata[ncdata.black_ind == 1].pb_75.mean()
print("False Negative")
ncdata[ncdata.black_ind == 1].pnb_75.mean()
print("True Negative")
ncdata[ncdata.black_ind == 0].pnb_75.mean()
print("Precision")
ncdata[ncdata.pb_75 == 1].black_ind.mean()
print("Recall")
ncdata[ncdata.black_ind == 1].pb_75.mean()
print("Accuracy")
len(ncdata[ncdata.black_ind == ncdata.pb_75])/len(ncdata)

print("\n90% Threshold")
print("False Positive")
ncdata[ncdata.black_ind == 0].pb_90.mean()
print("True Positive")
ncdata[ncdata.black_ind == 1].pb_90.mean()
print("False Negative")
ncdata[ncdata.black_ind == 1].pnb_90.mean()
print("True Negative")
ncdata[ncdata.black_ind == 0].pnb_90.mean()
print("Precision")
ncdata[ncdata.pb_90 == 1].black_ind.mean()
print("Recall")
ncdata[ncdata.black_ind == 1].pb_90.mean()
print("Accuracy")
len(ncdata[ncdata.black_ind == ncdata.pb_90])/len(ncdata)

print("\n\nRecalibrated Analysis \n")

print("Area under ROC curve")
roc_auc_score(ncdata['black_ind'], ncdata['bhat'])

print("\n50% Threshold")
print("False Positive")
ncdata[ncdata.black_ind == 0].pbhat_50.mean()
print("True Positive")
ncdata[ncdata.black_ind == 1].pbhat_50.mean()
print("False Negative")
ncdata[ncdata.black_ind == 1].pnbhat_50.mean()
print("True Negative")
ncdata[ncdata.black_ind == 0].pnbhat_50.mean()
print("Precision")
ncdata[ncdata.pbhat_50 == 1].black_ind.mean()
print("Recall")
ncdata[ncdata.black_ind == 1].pbhat_50.mean()
print("Accuracy")
len(ncdata[ncdata.black_ind == ncdata.pbhat_50])/len(ncdata)

print("\n75% Threshold")
print("False Positive")
ncdata[ncdata.black_ind == 0].pbhat_75.mean()
print("True Positive")
ncdata[ncdata.black_ind == 1].pbhat_75.mean()
print("False Negative")
ncdata[ncdata.black_ind == 1].pnbhat_75.mean()
print("True Negative")
ncdata[ncdata.black_ind == 0].pnbhat_75.mean()
print("Precision")
ncdata[ncdata.pbhat_75 == 1].black_ind.mean()
print("Recall")
ncdata[ncdata.black_ind == 1].pbhat_75.mean()
print("Accuracy")
len(ncdata[ncdata.black_ind == ncdata.pbhat_75])/len(ncdata)

print("\n90% Threshold")
print("False Positive")
ncdata[ncdata.black_ind == 0].pbhat_90.mean()
print("True Positive")
ncdata[ncdata.black_ind == 1].pbhat_90.mean()
print("False Negative")
ncdata[ncdata.black_ind == 1].pnbhat_90.mean()
print("True Negative")
ncdata[ncdata.black_ind == 0].pnbhat_90.mean()
print("Precision")
ncdata[ncdata.pbhat_90 == 1].black_ind.mean()
print("Recall")
ncdata[ncdata.black_ind == 1].pbhat_90.mean()
print("Accuracy")
len(ncdata[ncdata.black_ind == ncdata.pbhat_90])/len(ncdata)


### EITC
print("\n\nEITC Population\n")

print("Area under ROC curve")
roc_auc_score(eitc_ncdata['black_ind'], eitc_ncdata['predicted_prob_black'])

print("\n50% Threshold")
print("False Positive")
eitc_ncdata[eitc_ncdata.black_ind == 0].pb_50.mean()
print("True Positive")
eitc_ncdata[eitc_ncdata.black_ind == 1].pb_50.mean()
print("False Negative")
eitc_ncdata[eitc_ncdata.black_ind == 1].pnb_50.mean()
print("True Negative")
eitc_ncdata[eitc_ncdata.black_ind == 0].pnb_50.mean()
print("Precision")
eitc_ncdata[eitc_ncdata.pb_50 == 1].black_ind.mean()
print("Recall")
eitc_ncdata[eitc_ncdata.black_ind == 1].pb_50.mean()
print("Accuracy")
len(eitc_ncdata[eitc_ncdata.black_ind == eitc_ncdata.pb_50])/len(eitc_ncdata)

print("\n75% Threshold")
print("False Positive")
eitc_ncdata[eitc_ncdata.black_ind == 0].pb_75.mean()
print("True Positive")
eitc_ncdata[eitc_ncdata.black_ind == 1].pb_75.mean()
print("False Negative")
eitc_ncdata[eitc_ncdata.black_ind == 1].pnb_75.mean()
print("True Negative")
eitc_ncdata[eitc_ncdata.black_ind == 0].pnb_75.mean()
print("Precision")
eitc_ncdata[eitc_ncdata.pb_75 == 1].black_ind.mean()
print("Recall")
eitc_ncdata[eitc_ncdata.black_ind == 1].pb_75.mean()
print("Accuracy")
len(eitc_ncdata[eitc_ncdata.black_ind == eitc_ncdata.pb_75])/len(eitc_ncdata)

print("\n90% Threshold")
print("False Positive")
eitc_ncdata[eitc_ncdata.black_ind == 0].pb_90.mean()
print("True Positive")
eitc_ncdata[eitc_ncdata.black_ind == 1].pb_90.mean()
print("False Negative")
eitc_ncdata[eitc_ncdata.black_ind == 1].pnb_90.mean()
print("True Negative")
eitc_ncdata[eitc_ncdata.black_ind == 0].pnb_90.mean()
print("Precision")
eitc_ncdata[eitc_ncdata.pb_90 == 1].black_ind.mean()
print("Recall")
eitc_ncdata[eitc_ncdata.black_ind == 1].pb_90.mean()
print("Accuracy")
len(eitc_ncdata[eitc_ncdata.black_ind == eitc_ncdata.pb_90])/len(eitc_ncdata)

print("\n\nRecalibrated Analysis \n")

print("Area under ROC curve")
roc_auc_score(eitc_ncdata['black_ind'], eitc_ncdata['bhat'])

print("\n50% Threshold")
print("False Positive")
eitc_ncdata[eitc_ncdata.black_ind == 0].pbhat_50.mean()
print("True Positive")
eitc_ncdata[eitc_ncdata.black_ind == 1].pbhat_50.mean()
print("False Negative")
eitc_ncdata[eitc_ncdata.black_ind == 1].pnbhat_50.mean()
print("True Negative")
eitc_ncdata[eitc_ncdata.black_ind == 0].pnbhat_50.mean()
print("Precision")
eitc_ncdata[eitc_ncdata.pbhat_50 == 1].black_ind.mean()
print("Recall")
eitc_ncdata[eitc_ncdata.black_ind == 1].pbhat_50.mean()
print("Accuracy")
len(eitc_ncdata[eitc_ncdata.black_ind == eitc_ncdata.pbhat_50])/len(eitc_ncdata)

print("\n75% Threshold")
print("False Positive")
eitc_ncdata[eitc_ncdata.black_ind == 0].pbhat_75.mean()
print("True Positive")
eitc_ncdata[eitc_ncdata.black_ind == 1].pbhat_75.mean()
print("False Negative")
eitc_ncdata[eitc_ncdata.black_ind == 1].pnbhat_75.mean()
print("True Negative")
eitc_ncdata[eitc_ncdata.black_ind == 0].pnbhat_75.mean()
print("Precision")
eitc_ncdata[eitc_ncdata.pbhat_75 == 1].black_ind.mean()
print("Recall")
eitc_ncdata[eitc_ncdata.black_ind == 1].pbhat_75.mean()
print("Accuracy")
len(eitc_ncdata[eitc_ncdata.black_ind == eitc_ncdata.pbhat_75])/len(eitc_ncdata)

print("\n90% Threshold")
print("False Positive")
eitc_ncdata[eitc_ncdata.black_ind == 0].pbhat_90.mean()
print("True Positive")
eitc_ncdata[eitc_ncdata.black_ind == 1].pbhat_90.mean()
print("False Negative")
eitc_ncdata[eitc_ncdata.black_ind == 1].pnbhat_90.mean()
print("True Negative")
eitc_ncdata[eitc_ncdata.black_ind == 0].pnbhat_90.mean()
print("Precision")
eitc_ncdata[eitc_ncdata.pbhat_90 == 1].black_ind.mean()
print("Recall")
eitc_ncdata[eitc_ncdata.black_ind == 1].pbhat_90.mean()
print("Accuracy")
len(eitc_ncdata[eitc_ncdata.black_ind == eitc_ncdata.pbhat_90])/len(eitc_ncdata)
